#include <iostream>
#include <cmath>
#include "fraction.h"

using namespace std;

const int MAXN = 100;

Fraction matrix[MAXN][MAXN + 1]; // 矩阵
int n;                           // 矩阵行数，即未知量个数

/* 打印矩阵 */
void printMatrix()
{
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j <= n; j++)
        {
            cout << matrix[i][j] << " ";
        }
        cout << endl;
    }
}

/* 高斯消元 */
bool gaussianElimination()
{
    for (int i = 0; i < n; i++)
    {
        /* 找到列中绝对值最大的元素 */
        double max = abs(matrix[i][i].todouble());
        int maxrow = i;
        for (int j = i + 1; j < n; j++)
        {
            if (abs(double(matrix[j][i].todouble())) > max)
            {
                max = abs(matrix[j][i].todouble());
                maxrow = j;
            }
        }

        /* 如果该元素为0，则无法进行消元 */
        if (max == 0)
        {
            return false;
        }

        /* 交换行，将绝对值最大的元素移到主对角线上 */
        if (maxrow != i)
        {
            for (int j = i; j <= n; j++)
            {
                Fraction temp = matrix[i][j];
                matrix[i][j] = matrix[maxrow][j];
                matrix[maxrow][j] = temp;
            }
        }

        /* 消元 */
        for (int j = i + 1; j < n; j++)
        {
            Fraction factor = matrix[j][i] / matrix[i][i];
            for (int k = i; k <= n; k++)
            {
                matrix[j][k] -= factor * matrix[i][k];
            }
        }
    }

    /* 回带 */
    for (int i = n - 1; i >= 0; i--)
    {
        matrix[i][n] /= matrix[i][i];
        matrix[i][i] = 1;
        for (int j = i - 1; j >= 0; j--)
        {
            matrix[j][n] -= matrix[j][i] * matrix[i][n];
            matrix[j][i] = 0;
        }
    }

    return true;
}

int factorial(int _n)
{
    if (_n == 1 || _n == 0)
    {
        return 1;
    }
    else
    {
        return _n * factorial(_n - 1);
    }
}

void computeAdamsBash(int index)
{
    /* 输入未知量个数和矩阵 */
    n = index;
    int p = index;
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            Fraction f(pow(j, i), factorial(i));
            matrix[i][j] = f;
        }
        Fraction f1(pow(n, i + 1), factorial(i + 1)), f2(pow(n - 1, i + 1), factorial(i + 1));
        Fraction f = f1 - f2;
        matrix[i][n] = f;
    }

    /* 消元并输出结果 */
    if (gaussianElimination())
    {
        cout << "s = " << index << "\t" << "p = " <<p << "\t";
        for (int i = n - 1; i >= 0; i--)
        {
            cout << matrix[i][n] << "\t";
        }
        cout << endl;
    }
    else
    {
        cout << "No solution." << endl;
    }
}

void computeAdamsMoulton(int index)
{
    /* 输入未知量个数和矩阵 */
    int s = index;
    int p = s + 1;
    n = p;
    for (int i = 0; i < p; i++)
    {
        for (int j = 0; j < p; j++)
        {
            Fraction f(pow(j, i), factorial(i));
            matrix[i][j] = f;
        }
        Fraction f1(pow(s, i + 1), factorial(i + 1)), f2(pow(s - 1, i + 1), factorial(i + 1));
        Fraction f = f1 - f2;
        matrix[i][n] = f;
    }

    /* 消元并输出结果 */
    if (gaussianElimination())
    {
        cout << "s = " << index << "\t" << "p = " <<p << "\t";
        for (int i = n - 1; i >= 0; i--)
        {
            cout << matrix[i][n] << "\t";
        }
        cout << endl;
    }
    else
    {
        cout << "No solution." << endl;
    }
}

void computeBackward(int index)
{
    /* 输入未知量个数和矩阵 */
    int s = index, p;
    p = s;
    n = s + 1;
    for (int j = 0; j <= s - 1; j++)
    {
        Fraction f(1);
        matrix[0][j] = f;
    }
    matrix[0][s + 1] = -1;
    for (int i = 1; i <= p; i++)
    {
        for (int j = 0; j <= s - 1; j++)
        {
            Fraction f(pow(j, i), factorial(i));
            matrix[i][j] = f;
        }
        Fraction f3(-pow(s, i - 1), factorial(i - 1));
        matrix[i][s] = f3;
        Fraction f(-pow(s, i), factorial(i));
        matrix[i][n] = f;
    }

    /* 消元并输出结果 */
    if (gaussianElimination())
    {
        cout << "s = " << index << "\t" << "p = " <<p << "\t";
        cout << matrix[n - 1][n] << "\t";
        cout << 1 << "\t";
        for (int i = n - 2; i >= 0; i--)
        {
            cout << matrix[i][n] << "\t";
        }
        
        cout << endl;
    }
    else
    {
        cout << "No solution." << endl;
    }
}

int main()
{
    cout << "--------------------Adams-Bashforth formulas--------------------" << endl;
    for (int i = 1; i <= 5; i++)
        computeAdamsBash(i);
    cout << "--------------------Adams-Moulton formulas--------------------" << endl;
    for (int i = 1; i <= 4; i++)
        computeAdamsMoulton(i);
    cout << "--------------------Backward differentiation formulas--------------------" << endl;
    for (int i = 1; i <= 4; i++)
        computeBackward(i);

    return 0;
}
